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INTRODUCTION 


1  . 

With  the  advent  of  large  digital  computers,  numerical 
modeling  has  taken  an  important  place  in  atmospheric  science. 

<  Constructing  a  numerical  model  often  involves  two  steps.  The 

first  step  is  to  formulate  a  set  of  continuous  equations  which 
represent  the  dynamics  and  the  physics  of  interest  in  the 
atmosphere.  The  second  step  is  to  discretize  these  equations 
so  that  they  may  be  solved  on  a  computer.  ^  '  * 

Broadly  speaking,  the  discretization  methods  for  solving  the 
continuous  equations  may  be  grouped  as  Eulerian  or  Lagrangian 
(fully  Lagrangian  or  semi -Lagrangian )  methods.  The  Eulerian 
methods  may  be  further  classified  according  to  the  techniques 
employed  for  the  integration  in  time  and  the  approximations  used 
for  the  spatial  derivatives.  The  integration  in  time  can  be 
either  explicit  or  implicit  (fully  implicit  or  semi-implicit) 
while  the  spatial  discretization  may  themselves  be  grouped  under 
two  general  headings:  series  expansion  methods  and  finite  differ¬ 
ence  methods.  The  Galerkin  procedure  is  often  involved  in  the 
series  expansion  methods.  Some  commonly  used  series  expansion 
methods  are  finite  element,  spectral  and  collocation.  In  the 
finite  difference  method,  the  flux  form  is  often  preferred  to 
the  advective  form  in  order  to  maintain  in  the  discretized 
equations  integral  constraints  of  the  continuous  atmosphere. 

In  addition,  the  variables  are  staggered  around  the  grid  points 
»  for  the  sake  of  better  geostrophic  adjustment.  Of  the  staggered 

grids  B,  C,  D  and  E  (Arakawa  and  Lamb,  1977),  the  C-grid  has 

ft 

found  favor  amongst  numerical  modelers  (Arakawa  and  Lamb,  1977; 
Schoenstadt,  1978)  despite  criticisms  raised  by  some  authors 
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(Mesinger,  1973).  There  are  very  few  finite  difference  models 
that  employ  a  non-staggered  grid.  The  C-grid  is  currently  used 
in  the  Navy  Operational  Regional  Atmospheric  Prediction  System 
( Hodur ,  1982,  1987). 

There  are  many  problems  in  environmental  science  where  a 
central  concern  is  the  manner  in  which  a  trace  constituent  or 
water  vapor  is  transported  by  moving  fluid.  It  is  important  to 
find  a  suitable  finite  difference  algorithm  on  a  staggered  grid 
for  this  problem.  In  this  report,  we  use  a  simple  water  vapor 
transport  problem  (i.e.,  one-dimensional  linear  advection  of 
water  vapor)  to  study  the  advective  processes  on  the  staggered 
grid.  It  is  illustrated  that  caution  must  be  taken  in  defining 
flux  on  the  staggered  C-grid.  It  is  shown  that  accuracy  can  be 
lost  in  the  fourth-order  centered  differencing  if  flux  is  defined 
improperly.  A  consistent  way  of  defining  flux  with  the  fourth- 
order  finite  differencing  is  presented. 

2.  MODEL  PROBLEM  AND  DISCRETIZATION  SCHEMES 

To  illustrate  the  points  mentioned  above,  we  use  finite 
difference  schemes  to  discretize  the  one-dimensional,  constant 
coefficient  (u  =  const.)  advection  equation 
3q  3q 

—  +  u  —  =0  (2.1) 

9t  3x 

in  the  periodic  domain  [-1,1]  with  the  initial  condition 

2 

x 

q(x,t  =  0)  =  exp  l-  ( - )  ]  .  (2.2) 

0.2 

This  problem  is  the  simplest  prototype  of  a  model  involving  wave 
or  advective  processes.  For  the  test  purpose,  Eq .  (2.1)  can  be 
easily  written  in  flux  form  as 


2 


9q  OF 

—  +  —  =0 
3 1  3x 


(2.3a) 


F  =  uq  (2.3b) 

and  solved  on  the  staggered  grid  shown  in  Figure  1.  Note  that 
the  flux  is  defined  on  the  velocity  point.  Since  the  q  points 
and  flux  points  are  staggered,  we  shall  define  the  flux  in  terms 
of  surrounding  q  values  and  consider  the  spatial  derivative  by 
either  second-order  or  fourth-order  centered  differences. 

Schemes  (I)-(V),  which  involve  different  ways  to  define  the 
flux  and  to  approximate  the  spatial  derivative,  are  used  to 
discretize  (2.3).  The  definition  of  flux  and  space  differencing 
of  these  schemes  are  as  follows: 

Scheme  ( I ) 

3qi 

at 


=  -  6  vF 


x*  i 


(2.4a) 


<ui+l/2  +  '  ui +1/2 ' )  <ui+l/2  "  1 u i +1/2  1  > 

Fi  +1/2  =  -  <3i  +  -  <3 i  +1 


(2.4b) 


2.0 


2.0 


Scheme  (II) 


3qi  9  1 

-  -  (—  Fi  ”  —  ^ 3x  Fi ) 

3t  8  8 

(ui+l/2  +  1 UI +1/ 2 1 >  ^ ui +1/2  "  iui+l/2!) 

Fi +1/2  =  3i  +  3i+l' 


(2.5a) 


(2.5b) 


2 . 0 


2.0 


Scheme  (III) 
9<3i 


at 


-  -  <fj 


<Ji  +  3i+l 

Fi +1/2  =  -  ui +1/2 

2.0 


(2.6a) 


(2.6b) 


3 


F  F  F  F 


Figure  1.  Staggered  grid  of  two  variables. 


A 


Scheme  (IV) 


aqi  9  1 

-  =  -  (~  <fx  Fi  ~  ~  *3x  Fi)  (2.7a) 

dt  8  8 


<3i  +  Qi  +  1 

Fi +1/2  =  -  ui+l/2 

2.0 


(2.7b) 


Scheme  (V) 

3qi  9  1 

-  =  -  (-  <fx  Fi - <f3x  Fi)  (2.8a) 

3t  8  8 

Fi+l/2  =  l 9/16 (qi+i  +  qi)  -  l/16(qi+2  +  qi-iH  ^i+1/2  .  (2.8b) 

where 

Fi +1/2  -  Fi -1/2  Fi +3/2  -  Fi -3/2 

6X  Fi  -  and  tf3X  Fi  - 

AX  3  AX 

Forward  time  integration  is  used  in  Scheme  (I)  and  (II) 
while  the  second-order  leapfrog  time  integration  is  used  in 
Schemes  (III),  (IV)  and  (V)  for  their  time  discretization.  To 
suppress  the  computational  mode  in  the  leapfrog  time  integration, 
an  Asselin  filter  with  a  damping  coefficient  of  0.02  is  used. 

Note  that  the  first-order  upstream  flux  is  used  in  both  Schemes 
(I)  and  (II),  while  the  second-order  flux  is  used  in  Schemes 
(III)  and  (IV).  The  flux  in  Scheme  (V)  is  fourth-order  in 
accuracy  according  to  the  Taylor  expansion.  The  spatial 
derivatives  are  approximated  by  either  second-order  centered 
differencing  (Schemes  (I)  and  (III))  or  fourth-order  centered 
differencing  (Schemes  (II),  (IV)  and  (V)).  Thus,  Schemes  (III) 
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and  (V)  define  the  flux  and  approximate  the  3patial  derivative 
with  the  same  order  of  accuracy  (i.e.,  second-order  for  (III) 


and  fourth-order  for  (V)).  In  the  linear  advection  equation, 
space  differencing  of  Scheme  (I)  is  the  same  as  the  first-order 
upstream  differencing.  Note  that  Schemes  (II)  and  (IV)  involve 
a  mixture  of  accuracy  in  defining  flux  and  approximating  the 
spatial  derivatives.  The  fourth-order  Scheme  (V),  which  is 
for  the  C-grid  even  in  a  two-dimensional  situation,  is 
analogous  to  the  scheme  described  by  Gerrity  et  al.  (1972)  for 
the  non-staggered  grid  and  Campana's  (1973)  implementation  for 
the  B-grid. 

Recently,  Professor  A.  Arakawa  (personal  communication), 
in  the  context  of  solving  the  continuity  equation  in  a  primitive 
equation  model  using  isentropic  surfaces  as  a  vertical  coordinate, 
has  proposed  a  generalization  of  Takacs *  (1985)  scheme  which  has 
very  small  computational  dispersion  and  which  guarantees  positive 
results.  Arakawa's  scheme  is  third-order  in  accuracy.  We  will 
include  Arakawa's  scheme  in  our  results  for  comparison.  Formula¬ 
tion  of  Arakawa's  scheme  is  presented  in  Appendix  A.  In  addition 
to  Arakawa's  method,  we  will  also  test  the  positive  definite 
advection  scheme  of  Smolarkiwicz  (1983).  Similar  to  Arakawa's 
method,  Smolark iwicz ' s  method  involves  pred ictor -cor rector 
procedure;  but,  the  cost  of  Smolark iwicz ' s  method  is  lower 
compared  to  Arakawa's  method.  Detail  of  Smolark iwicz ' s  scheme 
is  in  Appendix  B. 
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MODEL  PROBLEM  ANALYSIS 


For  the  one-dimensional  constant  coefficient  advection  equa¬ 
tion,  Scheme  (I)  employs  the  first-order  upstream  differencing 
while  Scheme  (III)  employs  the  second-order  centered  differencing. 
Detailed  analysis  of  Scheme  (I)  and  Scheme  (III)  for  the  one¬ 
dimensional  linear  advection  equation  can  be  found  in  Haltiner  and 
Williams  (1980).  To  better  understand  the  accuracy  and  stability 
properties  of  Schemes  (IV)  and  (V)  we  now  follows  an  argument 
similar  to  that  given  by  Haltiner  and  Williams  (1980).  Consider 
the  following  fourth-order,  space-centered  differencing  and 
leapfrog  time  differencing  for  the  constant  coefficient  advection 
(Eq.  (2.1)). 


n+1  n-1 

9m  ~  “3m 


2  At 

n 

9m+  2 
+  b  - 


-  u  ( a 

9m- 2 


n  n 

9m+l  “  9m-l 


2ax 

n  n 

9m+3  "  9111-3 
+  c  - )  . 


(3.1) 


4 AX  6 AX 

Since  each  of  the  terms  with  coefficients  a,  b  and  c  are  valid 

approximations  for  3F/3x,  it  is  clear  that  (3.1)  will  be  a 

consistent  scheme  provided  that  a  +  b  +  c  =  1.  The  requirement 

for  (3.1)  to  be  fourth-order  accurate  is  that  the  second-order 

truncation  terms  vanish.  This  implies  that  the  sum  of  the 

Taylor  series  terms 

33q  3  93q  3 

a  ( — r  axj/3!(2ax),  b( — r-)  (2AX)  /3!(4AX), 

3x  dx 

and  a  similar  term. 


3  9  ■> 

c ( — 3)  (3ax)j/3! (6ax), 
3xJ 
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from  series  expansions  for  %4-l'  3m+2  and  ^m+3  ®ust  vanish. 

2  4 

Thus  the  resulting  truncation  error  in  (3.1)  is  0  ( At  )  +  0  (  ax  )  if 
a  +  4b  +  9C  =  0 .  (3.2) 

The  stability  will  be  investigated  next.  Using  the  usual 
notation,  ve  assume  a  solution  of  the  form 

n  .  ,  ianAt  +  iMmAx  ...  .. 

Qm  — AL  (3.3) 


and  define 

b  c 

f(MAx)  =  a  sin  pax  +  —  sin2  max  +  —  sin3  max  (3.4) 

2  3 

f’(MAx)  =  a  cos  Max  +  b  cos2  max  +  c  cos3  Max.  (3.5) 

Substitution  of  (3.3)  into  (3.1)  leads  to 

U  At 

sin  a At  =  -  -  f ( max )  .  (3.6) 

ax 

If  a  is  real,  neutral  solutions  result  with  no  damping  or  amplifi¬ 
cation,  as  evident  from  (3.3).  The  condition  for  a  real  is  that 
the  right  side  of  (3.6)  has  magnitude  less  than  or  equal  to  1, 
otherwise  it  is  complex.  To  insure  stability  for  all  wavelengths, 
it  is  necessary  to  find  the  maximum  magnitude  of  f ( max ) .  Let 
f((MAx)max)  be  the  maximum  value  of  f ( max ) .  Then  the  following 
criterion  must  be  satisfied  to  maintain  stability  for  all  wave 
numbers . 

u-t 

I - 1  <  1/If ( (MAx)maxI t  .  (3.7) 

AX 
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The  phase  speed  and  group  velocity  resulting  from  the  discretiza¬ 
tion  of  (3.1)  can  also  be  derived.  The  phase  velocity  Cp  is 

^  -1 

CF  =  u  *  -  sin  [  CR*  f  (max  )  ]  (3.8) 

MAX  CR 

while  the  group  velocity  Cg  is 

Cg  =  u  *  f • (MAx)/(l-[CR  *  f(MAx)]2)1/2  (3.9) 

where  CR  is  the  Courant  number. 

For  the  one-dimension  linear  advection  problem,  we  car 
easily  rearrange  the  leapfrog  time  discretization  of  (2.7)  and 
(2.8)  into  the  form  of  (3.1).  We  will  get  a  =  13/12,  b  =  -  1/12 
and  c  =  0  for  Scheme  (IV)  and  a  =  87/64,  b  =  -3/8  and  c  =  1/64 
for  Scheme  (V).  Note  that  for  the  non-staggered  grid  with 
fourth-order  difference,  as  given  in  Haltiner  and  Williams 
(1980),  a  =  4/3,  b  =  -1/3  and  c  =  0.  Although  Schemes  (IV)  and 
(V)  are  consistent.  Scheme  (IV)  is  not  a  fourth-order  scheme 
because  it  violates  (3.2).  Note  that  the  minimum  resolvable 
wavelength  of  2ax  is  stationary  in  both  Schemes  (IV)  and  (V) 
according  to  (3.8).  Also  there  are  physical  and  computational 
modes  in  the  numerical  solutions  as  revealed  by  (3.8). 

The  group  velocity  for  a  2ax  length  in  Scheme  (V)  is 
strongly  negative  with  Cg  =  -7/4u  while  the  group  velocity  for 
a  4ax  length  with  CR  =  0.7  is  almost  the  same  as  u.  It  can  be 
concluded  therefore  that  Scheme  (V)  can  lead  to  more  rapid 
spreading  of  noise  to  the  phase  velocity  in  the  very  short 
wavelengths  from  any  source  in  the  model,  physical  or  computa¬ 
tional.  Table  1  gives  the  ratio  of  finite  difference  wave  speed 
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Table  1. 

Cp/ U 

as  functions 

of  CR  and 

AX  . 

Schemes 

CR 

2  AX 

4  AX 

6  AX 

8  AX 

10  AX 

1 2D  Ax 

Non-staggered 

4th  Order 

a  =  4/3 

0.2 

0 

0.86 

0.97 

0.99 

1.00 

1.00 

b  =  -1/4 

0.4 

0 

0 . 89 

0.99 

1.00 

1.01 

1.01 

c  =  0 

0.6 

0 

0.98 

1.03 

1.03 

1.02 

1.01 

(IV) 

a  =  13/12 

0.2 

0 

0.69 

0.86 

0.92 

0.96 

0.98 

b  =  -1/12 

0.4 

0 

0.71 

0 .88 

0.93 

0.96 

0.98 

c  =  0 

0.6 

0 

0.75 

0.91 

0.95 

0.98 

0.99 

(V) 

a  =  87/64 

0.2 

0 

0.87 

0.97 

0.99 

0.99 

0.99 

b  =  -3/8 

0.4 

0 

0.91 

0.99 

1.00 

1.00 

1.00 

c  =  1/64 

0.6 

0 

1.00 

1.04 

1.03 

1.01 

1.01 

to  true  wave  speed,  Cp/u,  as  a  function  of  wavelength  in  terms 
of  ax  versus  Courant  number  CR  for  Schemes  (IV),  (V)  and  the 
non-staggered  fourth-order  finite  difference  for  the  constant 
coefficient  advection  equation.  From  Table  1,  we  conclude 
that  Scheme  (IV)  is  only  of  second-order  accuracy  despite  the 
additional  calculation  performed  for  the  fourth-order  spatial 
derivative  and  Scheme  (V)  is  a  proper  implementation  of  the 
fourth-order  finite  difference  scheme  on  the  staggered  grid 
shown  in  Figure  1. 

In  view  of  the  significant  accuracy  improvement  in  the 
fourth-order  schemes,  as  compared  to  the  second-order  schemes,  « 

we  now  concentrate  on  the  stability  analysis  of  the  fourth-order 
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Scheme  (V).  With  the  parameters  a,  b  and  c  in  Scheme  (V),  we 
found  that  (MAx)max  %  75°  and  the  stability  condition  (3.7)  of 
Scheme  (V)  gives 

UAt 

ICR  I  =  I - 1  i  0.71  .  (3.10) 

Ax 

Equation  (3.10)  indicates  that  the  CFL  criterion  associated  with 
Scheme  (V)  is  more  restrictive,  and,  for  a  given  u  and  ax,  the 
time  step  would  have  to  be  29  percent  smaller  for  the  fourth- 
order  Scheme  (V)  versus  the  second-order  Scheme  (III).  Note 
that  the  CFL  criterion  of  Scheme  (V)  is  the  same  as  the  non- 
staggered  fourth-order  finite  difference  of  Haltiner  and 
Williams  (1980) 

4.  NUMERICAL  RESULTS 

A  very  small  time  step  is  used  in  the  time  integrations 
presented  in  this  section  so  that  the  error  in  the  computation 
is  dominated  by  the  spatial  discretization  error  due  to  finite 
difference  methods. 

Figure  2  shows  the  approximate  solutions  at  t  =  2.0  obtained 
by  Schemes  (I)  -  (V)  and  the  scheme  of  Arakawa  with  number  of 
grid  points  N  =  32.  Note  that  Schemes  (I)  and  (II)  produces 
practically  identical  results  when  compared  with  Schemes  (III) 
and  (IV).  It  is  clear  from  Figure  2  that  Scheme  (V)  gives  a  much 
better  solution  than  do  Schemes  (III)  and  (IV).  Schemes  (I), 
(II),  and  the  Smolarkiwicz  scheme  and  Arakawa  scheme  generates 
positive  definite  fields.  The  Arakawa  scheme  is  associated  with 
the  least  damping  among  all  the  positive-definite  schemes. 
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Figure  2.  Analytical  and  numerical  solutions  with  N=32 
at  t=2 . 
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Figure  3  Illustrates  the  numerical  solutions  at  t  =  2.0 
obtained  by  all  the  positive-definitive  schemes.  Smolark ivicz ' s 
scheme  is  tested  with  various  S  factors  and  with  two  corrective 
steps.  There  seems  to  be  a  slightly  upwind  phase  shift  with  the 
Smolarkiwicz  scheme  for  S  =  1.  Moreover,  the  choice  of  S  factor 
does  not  seem  to  be  straightforward.  Both  the  Arakawa  and 
Smolarkiwicz  schemes  preserve  the  mass  very  well. 

Figure  4  shows  the  corresponding  L2  error  as  a  function  of 
the  number  of  grid  points  for  t  =  2.0. 

IN  1/2 

Here  L2  error  =  {—  E  lu(Xj,t)  -  ua(Xj,t)]  } 

N  j=l 

with  ua(Xj,t)  is  the  analytical  solution  and  u(Xj,t)  is  a 
computed  solution  at  grid  Xj  at  time  t.  The  different 
convergence  properties  of  different  schemes  can  be  easily  seen 
from  Figure  4.  For  example.  Scheme  (I)  and  (II)  converge  to  the 
analytical  solution  with  first-order  accuracy  (as  expected  from 
upstream  differencing)  while  Schemes  (III)  and  (IV)  converge  to 
the  analytical  solution  with  the  same  second-order  accuracy.  The 
fourth-order  accuracy  of  Scheme  (V)  is  obvious  when  compared  to 
Arakawa's  third-order  scheme.  Figure  3  reinforces  the  observa¬ 
tion  from  Figure  2  that  higher  accuracy  is  associated  with  the 
Scheme  (V)  solution.  In  particular,  the  Scheme  (V)  solution  is 
about  ten  times  more  accurate  than  the  Scheme  (III)  and  (IV) 
solutions  with  the  same  number  of  grid  points.  Because  Scheme 
(V)  only  requires  less  than  three  times  as  much  computation  as 
Scheme  (III),  there  is  a  net  gain  in  accuracy  by  using  fourth- 
order  Scheme  ( V ) . 
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Figure  3.  Analytical  and  numerical  solutions  of  positive 
definite  schemes  with  N*32  at  t=2. 
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L2  Error  at  t=2.</> 


5.  CONCLUDING  REMARKS 

In  this  report,  we  used  the  one-dimensional  constant 
advection  of  a  Gaussian  function  to  test  four  positive-definite 
schemes  ((I),  (II),  Smolarkivicz  scheme,  and  Arakawa  scheme)  and 
three  even-order  accurate  schemes  ((III),  (IV)  and  (V))  on  the 
staggered  grid.  It  is  shown  that  for  the  test  problem.  Schemes 
(I)  and  (II)  produce  identical  numerical  results.  The  Arakawa 
scheme  produces  the  least  damping  among  all  the  positive-definite 
schemes;  however,  there  are  expensive  computations  required  with 
the  Arakawa  scheme.  The  Smolarkiwicz  scheme  seems  to  be  cost 
effective.  The  slightly  upstream  shifted  phase  in  numerical 
solution  as  well  as  the  determination  of  the  S  factor  may  be 
the  drawback  of  the  scheme.  Of  the  even-order  accurate  schemes. 
Scheme  (V)  produces  very  accurate  results.  Stability  analysis 
indicates  that  a  reduction  in  time  step  size  (29%  smaller)  is 
required  when  Scheme  (V)  is  used  instead  of  Scheme  (III). 

Scheme  (V),  however,  produces  results  ten  times  more  accurate 
than  Schemes  (III)  and  (IV)  with  the  same  number  of  grid  points. 
Also,  Scheme  (V)  only  requires  less  than  three  times  as  much 
computation  as  Scheme  (III). 

Based  on  the  calculations  in  this  report,  we  suggest  that 
Schemes  (II)  and  (IV)  should  never  be  used  in  a  staggered  grid 
model.  This  is  because  the  accuracy  associated  with  the  flux  is 
inconsistent  with  the  accuracy  of  the  approximations  to  the 
spatial  derivatives.  The  additional  calculations  performed  in 
the  fourth-order  centered  differencing  in  these  two  schemes  do 
not  increase  model  accuracy. 
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Finally  we  recommend  that  for  the  purpose  of  horizontal 
advection  of  water  vapor  mixing  ratio  on  the  C-grid,  Scheme  (V) 
should  replace  (IV)  in  the  Navy  Operational  Regional  Atmospheric 
Prediction  System  (Hodur  1982,  1987).  Improvement  of  the  model 
accuracy  from  the  replacement  is  expected. 

In  closing  we  note  that  there  are  many  other  algorithms  that 
can  be  applied  to  the  transport  of  trace  constituent  or  water 
vapor  in  atmospherical  models.  A  more  complete  review  can  be 
found  in  Rood  (1987).  In  addition  to  many  useful  finite 
difference  schemes,  there  are  Lagrangian  and  series  expansion 
Eulerian  methods.  The  pseudospectral  (collocation)  and  the 
semi-Lagrangian  methods  are  of  special  interest.  Rasch  and 
Williamson  (1989)  used  shape  preserving  interpolation  schemes 
in  the  semi-Lagrangian  method  to  maintain  the  local  monotonicity 
of  the  simulated  fields.  Since  the  semi-Lagrangian  method  can  be 
used  in  conjunction  with  the  finite  difference  -  finite  element 
and  spectral  approaches  with  a  larger  time  step  than  the  Eulerian 
methods  -  the  shape-preserving  semi-Lagrangian  scheme  seems  to  be 
very  useful . 

The  efficiency  of  the  spectral  or  collocation  method 
often  lies  in  its  exponential  convergence  property  (Fulton 
and  Schubert,  1987).  For  the  linear  advection  equation,  the 
pseudospectral  method  is  the  same  as  the  Galerkin  or  tau 
method  except  that  the  2&x  wave  is  stationary  in  the  pseudo¬ 
spectral  method.  Finite  differences  on  non-staggered  grids  and 
Chebyshev  collocation  methods  have  been  applied  to  the  same 
advection  problem  studied  in  this  report  (Fulton  and  Schubert 
(1987)  )  . 
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APPENDIX  A 


Arakawa's  method  is  in  the  prediction-correction  form. 
Predictor: 

-  At  <fxFj.  (Al) 

where 


Fi  +  l/2  =  u  i  +  1/2  <3i  +  ui  +  l/2  q?  +  l 


( A2 ) 


ui +1/2  =  °-5  ( ui  +1/2  +  lui+l/2^ 


(A3) 


ui  +  1/2  =  °-5  ( ui +1/ 2  ‘  *  u i  +  1/2  1  > 


( A4  ) 


Corrector : 


n+1  n  . .  .  * 

q i  =  Qi  “  <*x  Fi 


( A5 ) 


where 
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ui+l/2 


*  n 
qj  +  1  +  <3i 
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+  ui  +1/2 


*  ,  n 

qj  *  qj+i 

2.0 
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Gi+l/2  =  “  “i+1/2  tui+l/2  ^ i+1/2  Ui  +  1  "  q?J  ~ 


~  ~  +  * 
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-  2qi  +  qi+i) 

(All) 

r i  +  1/2  ~ 
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Here  the  superscript  n  indicates  the  nth  time  level  during  time 

-15 

integration  and  e  =  10 
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APPENDIX  B 


The  Smolarkivicz  scheme  is  in  the  pred ictor -corrector  form. 
Predictor: 


(ui+l/2  +  lui+l/2l)  n  (ui+l/2  "  lui+l/2|)  n 

(Bl) 

Fi+l/2  = 

4i  x 

2.0  2.0 

4l+l 

* 

<3i  = 

Qi  -  At  <fx  F i 

( B2  ) 

Corrector 

; 

Ax  2  *  * 

(lui+1/2l -  At  ui+1/2)(qi+i  -  Qi) 

2.0 

_  *  c 

( B3  ) 

u  i  +  1/2  = 

Ax  *  * 

- (Qi  +  Qi+1  +  «) 

2.0 

(B4) 

* 

<  ui  +  1/2  +  1 ui +1/2 1 )  *  ( ui +1/2  * 

1 u i +1/2  1  ) 

* 

Fi+l/2  = 

<3i  4 

2.0 

2.0 

'll  +1 

n  +  1 

<3i 

=  q i  -  At  <fx  Fj 

( B5 ) 

Here  <=  =  10  S  is  an  engineer  factor  used  to  improve  the 

quality  of  solution  from  experiment.  The  superscript  n  indicates 
the  nth  time  level  during  time  integration.  The  corrector  step 
can  be  repeated  to  improve  the  solution  according  to  Smolarkivicz 
(1983).  We  have  also  tested  the  scheme  with  two  corrective 
s  teps . 
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